##############################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)

##############################################################################

option_list <- list(
	make_option(c("--class_type"), type = "character") ,
  make_option(c("--input_file"), type = "character") ,
  make_option(c("--images_path"), type = "character")
)

if(1!=1){
  
  class_type <- "DGC"
  work_dir <- "~/20220915_gastric_multiple/dna_combinePublic"
	input_file <- paste0( work_dir , "/images/Con_Exclusive/MutuallyExclusive.DGC.tsv")
	images_path <- paste(work_dir,"/images/Con_Exclusive/",sep="")
}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

input_file <- opt$input_file
images_path <- opt$images_path
class_type <- opt$class_type

###########################################################################################

dat_rtk <- data.frame(fread(input_file))

###########################################################################################

gene_list <- subset( dat_rtk , p < 0.05 )
gene_list <- unique(c( gene_list$geneA , gene_list$geneB ))
dat_rtk <- subset( dat_rtk , geneA %in% gene_list & geneB %in% gene_list )

###########################################################################################

dat <- dat_rtk
dat$OR2 <- ifelse( dat$OR > 1 , dat$OR , -1/dat$OR )
dat$OR2 <- ifelse( dat$OR2 > 4 , 4 , dat$OR2 )
dat$OR2 <- ifelse( dat$OR2 < -4 , -4 , dat$OR2 )
dat$OR2 <- ifelse( dat$p > 0.05 , 0 , dat$OR2 )

#dat$p <- ifelse( dat$p < 0.001 , 0.001 , dat$p )
#dat$p2 <- -log10(dat$p)
#dat$p2 <- ifelse( dat$OR > 1 , dat$p2 , -dat$p2 )

###########################################################################################
## 每个基因都有对应值
## 产生矩阵
gene_list <- unique(c(dat$geneA , dat$geneB))

dat_matrix_p <- matrix("" , ncol = length(gene_list) , nrow = length(gene_list) )
colnames(dat_matrix_p) <- gene_list
rownames(dat_matrix_p) <- gene_list

dat_matrix_or <- matrix("" , ncol = length(gene_list) , nrow = length(gene_list) )
colnames(dat_matrix_or) <- gene_list
rownames(dat_matrix_or) <- gene_list

for(i in 1:length(gene_list)){
	for( j in 1:length(gene_list)){
		if(i!=j){
			or <- subset( dat , 
				(geneA == gene_list[i] & geneB == gene_list[j] ) |
				(geneB == gene_list[i] & geneA == gene_list[j] )
				)$OR2
			p <- subset( dat , 
				(geneA == gene_list[i] & geneB == gene_list[j] ) |
				(geneB == gene_list[i] & geneA == gene_list[j] )
				)$p
			if(length(or)!=0){
				dat_matrix_or[i,j] <- as.numeric(or)
				dat_matrix_p[i,j] <- as.numeric(p)
			}else{
				dat_matrix_or[i,j] <- 0
				dat_matrix_p[i,j] <- 1
			}
		}else{
				dat_matrix_or[i,j] <- 0
				dat_matrix_p[i,j] <- 1
			}
	}
}

dat_matrix_or2 <- apply( dat_matrix_or , 1 , as.numeric)
rownames(dat_matrix_or2) <- gene_list
colnames(dat_matrix_or2) <- gene_list

dat_matrix_p2 <- apply( dat_matrix_p , 1 , as.numeric)
rownames(dat_matrix_p2) <- gene_list
colnames(dat_matrix_p2) <- gene_list

get_upper_tri <- function(cormat){
  cormat[lower.tri(cormat)]<- NA
  return(cormat)
}
upper_tri_or <- get_upper_tri(dat_matrix_or2)
melted_cormat_or <- melt(upper_tri_or , na.rm = TRUE)

upper_tri_p <- get_upper_tri(dat_matrix_p2)
melted_cormat_p <- melt(upper_tri_p , na.rm = TRUE)
colnames(melted_cormat_p)[3] <- "p"

###########################################################################################
## 
melted_cormat_or$id <- paste0( melted_cormat_or$Var1 , "_" , melted_cormat_or$Var2 )
melted_cormat_p$id <- paste0( melted_cormat_p$Var1 , "_" , melted_cormat_p$Var2 )

result_use <- merge( melted_cormat_or , melted_cormat_p[,c("id" , "p" )] , by = "id" )
result_use$p_text=ifelse(result_use$p>=0.05,"","*")
result_use$p_text=ifelse(result_use$p<0.05 & result_use$p>0.01,"*",result_use$p_text)
result_use$p_text=ifelse(result_use$p<0.01 & result_use$p>0.001,"**",result_use$p_text)
result_use$p_text=ifelse(result_use$p<0.001 ,"***",result_use$p_text)

result_use$Var1 <- as.character(result_use$Var1)
result_use$Var2 <- as.character(result_use$Var2)
result_use$Var1 <- factor(result_use$Var1 , levels = gene_list , order = T )
result_use$Var2 <- factor(result_use$Var2 , levels = gene_list[length(gene_list):1] , order = T )

p <- ggplot(data = result_use, aes(Var1, Var2, fill = value))+
  geom_tile(color = "white")+
  scale_fill_gradient2(
  	low = rgb(60,83,164 , maxColorValue = 255), 
  	high = rgb(213,73,65 , maxColorValue = 255), 
  	mid = rgb(251,251,251 , maxColorValue = 255), midpoint = 0 , space = "Lab", 
  	name="Odd Ratio" , breaks = c(-4,-2,0,2,4) , labels = c(">=4(Mututally exclusive)" , "2" , 0 , 2 , ">=4(Co-occurence)") ) +
  scale_x_discrete(position = "top") +
  geom_text(aes(x = Var1 , y = Var2 , label=p_text),size= 5 , color = "black") + 
  theme_minimal()+ 
  theme(axis.text.x = element_text(angle = 90, vjust = 1, 
                                   size = 12, hjust = 0))+
  scale_y_discrete(limits=rev) + ## y轴翻转
  coord_fixed() +
  theme(
    legend.position = c(0.8,0.2),
    legend.title = element_text(size = 12,color="black",face='bold'),
    panel.grid.major=element_blank(),
    panel.grid.minor=element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(),
    legend.text = element_text(size = 8,color="black",face='bold'),
    axis.text.y = element_text(size = 16,color="black",face='bold'),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(size = 16,color="black",face='bold') ,
    axis.line = element_line(size = 0.5)) 

out_name <- paste0( images_path , "/conOccurrence_exclussive.",class_type,".pdf" )
ggsave( out_name , p )